Synchronization transition of heterogeneously coupled oscillators 

on scale-free networks 
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We investigate the synchronization transition of the modified Kuramoto model where the oscilla- 
tors form a scale-free network with degree exponent A. An oscillator of degree fc, is coupled to its 
neighboring oscillators with asymmetric and degree-dependent coupling in the form of Jk^~ . By 
invoking the mean-field approach, we determine the synchronization transition point J c , which is 
zero (finite) when 77 > A — 2 (77 < A — 2). We find eight different synchronization transition behaviors 
depending on the values of 77 and A, and derive the critical exponents associated with the order pa- 
rameter and the finite-size scaling in each case. The synchronization transition is also studied from 
the perspective of cluster formation of synchronized vertices. The cluster-size distribution and the 
largest cluster size as a function of the system size are derived for each case using the generating 
function technique. Our analytic results are confirmed by numerical simulations. 
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I. INTRODUCTION 

Synchronization of oscillations is one of the fundamen- 
tal nonlinear phenomena in biology, physics, chemistry, 
communication science, and many other branches of sci- 
ence and engineering [lj . Recently, the dynamics of syn- 
chronization of oscillators located at each vertex in com- 
plex networks has attracted much attentions. That is 
because the small-world feature of complex networks is 
closely related to their synchronizability. By the small- 
world feature, we mean that the average separation (d) 
between a pair of vertices scales at most (d) ~ logiV, 
where N is the number of vertices in the system. It was 
shown that a small-world network model introduced by 
Watts and Strogatz Q is more synchronizable than reg- 
ular lattice However, such feature is not observed in 
scale-free (SF) networks. SF networks are the networks 
that exhibit a power-law degree distribution Pd(k) ~ k~ x 
and degree k is the number of edges connected to a given 
vertex In SF networks, the heterogeneity in the de- 
gree distribution suppresses their synchronizability ||. 
Thus, it was desired to introduce a new dynamic model 
that prompts SF networks to be more synchronizable. 

The dynamics of synchronization is described by vari- 
ous forms of coupled equations. A linearly coupled model 
is probably the simplest one. In the model, N oscil- 
lators are coupled when they are connected via edges. 
Coupling constant is normally symmetric; however, it is 
not necessarily symmetric to achieve a better synchro- 
nizability. This case can happen in SF networks: It was 
shown recently [|| that the synchronizability becomes 
maximum when information flow diffuses and reaches a 
uniform stationary state over the entire system. Here, the 
mapping from synchronization dynamics to information 
flow can be naturally introduced, because the linearly 
coupled equation is nothing but the diffusion equation. 
It was shown that the uniform-stationary state can be 



reached by introducing asymmetric and weighted cou- 
pling strength between a pair of vertices or oscillators. 

To be specific, the dynamic model with the asymmetric 
coupling strength is written as 
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for i = 1, . . . , N. Here, fa is the phase of an oscillator 
located at vertex i, f(<j>) describes the dynamics of an in- 
dividual oscillator, J is the overall coupling strength, fej 
is the degree of vertex i, and a^- is an element of the ad- 
jacent matrix, which is 1 if vertices i and j are connected 
and otherwise. h{4>i) is the output function and take a 
form of h(<pi) = 4>i f° r the linear case. It is noteworthy 
that the coupling strength of Eq. is asymmetric and 
weighted due to the factor l/k]~ v unless 77 = 1. When 
77 > 0, vertices with large degree can influence other ver- 
tices significantly on regulating phases due to their large 
numbers of connections; on the other hand, when 77 < 0, 
the influence is reduced. It was found [|| that the system 
is most synchronizable when 77 = 0, irrespective of the 
value of the degree exponent of a given SF network. 

In this paper, we study the pattern of synchroniza- 
tion transition for a modified Kuramoto model Q , which 
is similarly modified with the asymmetric and weighted 
coupling strength as 



d<t>i 
dt 



J 



N 

E 



di,- sin 



(2) 



The oscillators are located at each vertex i = 1,...,N 
of a SF network with degree exponent A. Here, Wj is 
the natural frequency of the ith oscillator selected from 
the Gaussian distribution g(u) = e~ w * 'I 2 j '\/2tt '. We find 
that the modified Kuramoto dynamic model displays a 
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very complex and rich behaviors in the space of the two 
tunable parameters (77, A). 

The synchronization transition from a desynchronized 
to a synchronized state occurs at the critical point J c . 
For small J -C J c , the coupling strength is so weak that 
an individual vertex maintains its own phase different 
from others; therefore, the entire system is desynchro- 
nized. As the coupling strength J increases, a cluster of 
vertices is more likely to be coupled, to be in a common 
or almost the same phase, and thus forms a cluster of 
synchrony. Size of such clusters becomes diverse as the 
coupling strength J increases. At the critical point J c , 
the system reaches a self-organized state, and the cluster- 
size distribution follows a power law, 



n{s) 



(3) 



in Sec. IIIII The size distribution of synchronized clusters 
and the largest cluster size at the critical point are solved 
in Sees. IIVI and El respectively. The finite-size scaling 
analysis for the order parameter is performed and the re- 
sults are checked numerically in Sec. IVII Summary and 
discussion follow in Sec. IVIII 



II. ORDER PARAMETER EQUATION 

In this section, we analyze the modified Kuramoto 
equation J5J in the framework of the mean-field approach 
by constructing a self-consistent equation for a local field. 
To proceed, we define fj and 6i as the amplitude and 
phase of the local field at vertex i, respectively, via 



in the thermodynamic limit. For J » J c , the power- 
law behavior no longer holds and the entire system is 
synchronized. 

The order parameter of the synchronization transition 
is defined as 



1 



re™ =-Ve*. 
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(4) 



In the synchronized state, the phases fa of each vertex 
are narrowly distributed around an average phase 9. The 
amplitude r of the order parameter has a finite value; on 
the other hand, r ~ in the desynchronized state. Thus, 
the exponent (3 associated with the order parameter is 
defined via the relation, 



(5) 



where A = ( J — J c )/J c . In finite-size systems, the order 
parameter is described in terms of a scaling function as 



r ~ iV-^V(AiV 1/ft ). 



(6) 



In the recent works (3, li| , the nature of the transitions 
and the finite-size scalings have been studied for the case 
of 77 = 1 . In this work, we determine the order parameter 
and the size distribution of synchronized clusters for gen- 
eral rj using the mean-field approach and the generating 
function technique. Moreover, we construct a finite-size 
scaling function for the order parameter, and determine 
the exponent /x. Even for a simple extension of r/ =/= 1, 
we find that the obtained result is very rich. There ex- 
ist eight distinct transition behaviors depending on the 
values of rj and A. Therefore, the result can be helpful 
in understanding diverse dynamic phenomena arising on 
SF networks. 

The paper is organized as follows: In Sec. ^ we first 
introduce and apply the mean-field approach to the dy- 
namic equation (0). We construct a self-consistent equa- 
tion for a local field and determine the order parameter. 
Next, the critical point is determined and the behavior 
of the order parameter near the critical point is obtained 



TV 



(7) 



Then, Eq. J2J) is rewritten in terms of the local field as 



— 7— = — JviW sin( 
dt y 



i). 



(8) 



Once the amplitude fj and the phase 6{ of the local field 
are determined, one can solve Eq. (JSJ easily. The local 
field fj is determined in a self-consistent manner. 

We consider the probability density pf\fauj)d(j) that 
the phase of an oscillator i with natural frequency lo lies 
between </> and 4> + d<fi in the steady state [TFj . Using a 

previous result [Toj that pf (fau))d(j) is inversely propor- 
tional to the speed of 4>, one can obtain that 



> 2-k\oj — 10* j sin(0 — 6i)\ 



if \oj\ < oj* 



otherwise, 



(9) 

where w*^ = Jr^fc^. This result implies that an oscillator 
i with natural frequency u> has its phase that is locked 
at <f> — 0~i + sin _1 (tj/(jj >l »,i) and dfa/dt = if \u>\ < 0J* t i. 
Otherwise, its phase drifts with a finite speed, dfa/dt 7^ 
0. Next, we can evaluate the order parameter using the 
stationary probability density in Eq. (JSJ as 



1 E 



dLug(u) I dfaD^ (faujy 4, . 



(10) 



Although fi, 9i, and pf\<f>\uj) can fluctuate over i in 
the steady state, we assume here that they depend only 
on degree ki. This is a mean field approximation. Keep- 
ing only the degree-dependent fluctuations, one can ob- 
tain a self-consistent equation for the local field through 
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Eq. 



as 



-~{k)e lS{k) = J2 P{k'\k) 



fc'=i 

/OO />27T 
dujg(uj) / #/j (s) (0|w,A;')e # , (11) 
-oo J0 

where (<f>\u, k) is given by the right hand side of 
Eq. © with oj*(k) = Jr(k)k ri replacing P(k'\k) 
denotes the probability that a neighboring vertex of a 
given vertex with degree k has degree k' , and k m is the 
natural cutoff of degree. Here, we consider only the case 
that the network is random and does not have any type 
of degree-degree correlation; then, P(k'\k) can be writ- 
ten as k'P d (k')/(k) with (k) = E fc fcP d (fc). After that, 
one can see that both f(k) and 9{k) are independent of 
degree fc,_and therefore, we can drop the fc-dependence 
in f and 9 from now on. 

The last integral of Eq. ltTT|l is evaluated as 



;,k) 



,i<t> 



i(w/w*(k)) — i\J (w/w»(fc)) 2 — 1, (u > w*(k)), 
< + V 1 ~ (w/w,(fc)) a , (H<w,(ft)), 

i(bj/uj lf (k)) + iyj (cj/a>»(fc)) 2 — 1, (w < — a>*(fc)). 

(12) 

The remaining integration in Ea. (|ll|l for uj > co t (k) and 
lj < LJ*(k) cancel out due to the fact g(co) = g(—uS). As 
a result, only the oscillators having the frequency within 
the range \u>\ < contribute to the local field in 

Eq. J2J. Thus, one obtains 



fr[ ( k ) V V^*( fc ) 



,(fc) 



U! 



(13) 



which is the self-consistent equation for f. Note that 
f is contained in = Jfk v . After the local field 

is obtained, the order parameter in Eqs. J3J or (|1(J|) is 
calculated as 



= Y J Pd{k) dug{d)K\ 

•/-w.(fc) v 



w*(/c) 



(14) 



= e-^'l 2 /^ to Eq. JT3J|, one can derive the local 
field f as 



2^ 7I\ ^*( fc ) 



fc=i 



*/ 2 . 1 
d<fi cos i 

tt/2 V27T 



[w«(fc) sin 0] 2 /2 



g (n-l/2)!(-l)"g^"- ! "'''- 1 (jf)2n+1 



i! (n + l)!2™+ 3 / 2 # 



A-l 



£i n (jf) 2 " +i , 



(15) 



n=0 



where we used the Taylor expansion of e -( UJ *( k ) sln <t>) / 2 
and the integration f^^ 2 d(f> cos 2 (f>sm 2n <fr = it l l 2 (n — 
1/2)!/ (2(n + 1)!). Similarly, the order parameter is eval- 
uated as 

~ (n- 1/2)! (-l)™^-"- 2 "" 



2 n+3/2 n | ( n+ l)| ffA 



2n+l 



(16) 



n=0 



If rj < 0, the generalized harmonic numbers in A n and 
A n are finite, and then they can be represented in terms 
of the Ricmann zeta functions for all n below and they 
are denoted as B n and B n , respectively. That is, 



-An 

and 



Br, 



A, 



B„, — 



(n - l/2)!(-l)"C(A - r) - 2m/ - 1) 
2"+ 3 / 2 n!(n+l)!C(A-l) 

(n - l/2)!(-l)"C(A - i] - 2m/) 



2«+3/2n!(n+l)!C(A) 



(17) 



(18) 



respectively Using these formulae, the local field and the 
order parameter are determined by Eqs. (|15[) and ((16(1 . 

On the other hand, it is remarkable that when < q < 
1, the generalized harmonic number diverges as H q n ~ 
(m+ l) 1 ^ q /(l-q) + ((q) + 0(m- q ) in the m -> oo limit, 
which is shown in Appendix. Then, Eqs. I|15|l and (|16|) 
are divided into analytic and singular parts as 

f = B n {Jr) 2n+1 + C(Jffc^)(Jf) (A - 2 )/" (19) 



and 



r = £ B n (Jf) 2n+1 + C(Jfkl)(Jf) 



(A-l)/r, 



(20) 



III. SYNCHRONIZATION TRANSITION 

In this section, we solve the self-consistent equation, 
Eq. HI 3(1 explicitly, and then investigate the behavior of 
the order parameter near the critical point via Eq. ((14(1 . 
To proceed, we first recall that the degree distribution is 
given in a closed form as Pd(k) = k~ x /Bt£ for A > 2, 
where Bt£ is the generalized harmonic number, defined 
by = YJk=i k ~ q and k m ^ TVV(A-i). Substituting 



respectively, where the functions C{x) and C(x) are de- 
fined in Appendix. In the x — > oo limit corresponding to 
the thermodynamic limit, C{x) and C(x) reduce to Coo 
and Coo, respectively, defined as 



Cry-, 



Coo — 



[(A-2t ? -2)/2t ? ]![(2-A-> ? )/27 ? ]! 

7?2 (A+4r,-2)/2r,[( A + ^ _ 2)/2t?]!C(A - 1) ' 

[(A-2 ?/ -l)/2 ?? ]![(l-A-r/)/2r/]! 
T?2 (A+^-i)/2r 3 [(A + 7/ - 1)/2t/]!C(A) ■ 



(21) 
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Thus, the local field and the order parameter are written 
as 

oo 

f = B n (Jf) 2n+1 + CUJr) {X ~ 2)h + ■■■, (22) 



and 



r = J2 B n (Jf) 2n+1 + C^Jf)^- 1 )/" + . .. , (23) 

n=0 

for Jfk^ 1. We remark that the singular terms appear 
only in the limit Jrfc^ — > oo. For the case of JffcJJ, <C 1, 
however, Eqs. I|15|l and (|16fl arc valid. 

Next, we determine the critical point. To proceed, we 
investigate the behavior of the local field as a function of 
J, which depends on the sign of ?y. 

(i) In the case of r] < 0, A n and A n are finite. 
One can see from Eq. I|15(l that the local field is zero 
for AqJ < 1 and non-zero for A$J > 1. The order pa- 
rameter behaves in the same manner as that of the local 
field from Eq. I|16[) . Thus, we obtain the critical point as 



Jc = T = 



1 2V2 H^- 1 



(24) 



As A — > 00, the critical point J c approaches 2y/2/y/n ~ 
1.60 in the limit N — > 00, which is consistent with that 
found in case of the globally-coupled oscillators 0. 

When J > J c , the local field f and the order parameter 
r are non-zero. When J is close to J r , 



and 



(25) 



(26) 



where A = ( J — J c )/J c . Thus, we obtain that (3 = 1/2. 
Again, this result is consistent with the one obtained from 
the globally-coupled oscillators 0. 

(ii) In the case of 77 > 0, the singular terms in 
Eqs. H22f) and (|23[l can be crucial in determining the criti- 
cal point and the order parameter. Depending on relative 
magnitude of A and 77, we divide the case of 77 > into 
four subcases: 

(I) When < 77 < (A - 2)/3 (i.e., A > 3?7 + 2), f ~ 
B Jf + BxJ 3 f 3 + ■ ■ ■ for small f from Eq. 
Then J c and f behave as those for 77 < presented 
in Eqs. Jg3J and 

(II) When (A - 2)/3 < r\ < A - 2 (i.e., 77 + 2 < A < 
377 + 2), the dominant contribution is made from 
the singular term of Eq. l(2"2|> . Then 



leading to 

Jc - l/Bo 

and 



2V2 C(A-l) 
0FC(A-7/-l) 



r ~ A*7/(*-2-i7)_ 



(28) 



(29) 



(III) When A-2<?7<A-1 (i.e., -q + 1 < A < 77 + 2), 
the critical point in Eq. (|24|l for finite N behaves 
as 



J ~ fc -(r;-A+2) _ Ar _(^_A+2)/(A-l)_ 



(30) 



Thus, it approaches zero in the thermodynamic 
limit, f is always positive unless J is zero as 
f ^ j(A-2)/(t7-A+2) for gmall J an d r - Jf - 

jr,/(r,-A+2)^ 

(IV) When 77 > A - 1 (i.e., A < 77 + 1), we obtain that 
r ~ ( Jf)^ -1 )/' 7 . Using the result of f obtained in 
(III), wc obtain that 



r ~ j(A-l)/(r;-A+2)^ 



(31) 



We summarize the result as follows: When 77 < A — 2 
(in the (I) and (II) cases), the critical point J c is finite; 
however, when 77 > A — 2 (in the (III) and (IV) cases), 
J c = in the thermodynamic limit N — ► 00. Thus, the 
critical exponent (3 associated with the order parameter 
is defined through the relation, r ~ A 13 (r ~ J 13 ) for the 
former (latter) case. The exponent (3 is evaluated in each 
case as follows: 



( 1/2 in (I), 

It7/(A-2-77) in (II), 

I 77/(7/ - A + 2) in (III), 

l(A-l)/fa-A + 2) in (IV). 



(32) 



f ~ B Q Jf +C 00 {Jf) 



(A-2)/r, 



(27) 



The result of the critical point is consistent with those 
of other phase transition problems such as the percolation 
transition and the epidemic spreading in SF networks. 
Moreover, the result for the case of 77 = 1 reduces to the 
previous result H,E|. Moreover, the result f3 = 1/2 for 
7/ = 1 and A > 5 is reduced to the mean-field result in 
regular lattice. 



IV. CLUSTER FORMATION OF 
SYNCHRONIZED OSCILLATORS 

In this section, we investigate in detail how the coupled 
oscillator system develops its synchrony as the coupling 
strength increases. To this end, we study the formation of 
clusters comprising synchronized vertices as a function of 
the coupling strength J. We use the generating function 
approach to derive the cluster-size distribution. 
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A. Cooperative versus background synchrony 

The order parameter averaged over the natural fre- 
quency distribution g(u>) can be written as 

n < E( c ° s2 & + sin2 ^) + s - &)) 



(33) 

from Eq.@. Here, the brackets represent the average 
over g(u). For the case of J = 0, each element oscillates 
independently, so that (cos(</>i —<j>j)) = f° r * 7^ J- Thus, 
the order parameter is evaluated as 



rj=o 



N 



(34) 



As J increases, clusters comprising synchronized oscilla- 
tors are more likely to form. We here define a cluster as a 
group of vertices (or oscillators) which arc connected and 
in the same coherent state: Two oscillators are coherent 
if its time-average correlation function C'ij , defined as 



(*1 



VJ t=t Q +l 



(cos(&(t)-^ (*))), (35) 



is larger than a preassigned threshold value Cth- We 
choose C t h value to generate the cluster-size distribu- 
tion in a power law form at the critical point J c . As 
such clusters form, the term of X)i S :J ^i( C0S ( < ^i — 4>j)) 
becomes nonzero and dominant. The order parameter is 
then evaluated as 



N 



(36) 



where k is the index of cluster and s K is the size of cluster 
K, i.e., the number of vertices within the cluster k. Note 
that J2 K s K = N, and Eq. JjESJ reduces to Eq. $Q when 
J = because each cluster size is 1. In the synchronized 
state when J ^> J c , the size of the largest cluster, de- 
noted as S, is of O(N), and thus the order parameter is 
approximately given as 



S/N. 



(37) 



Next, we study the cluster-size distribution and the size 
of the largest cluster as a function of J. 

The dynamics of cluster merging with increasing J 
results in the change of the cluster-size distribution. 
Let n(s) be the number of s-size clusters. Then 
y^ s sn(s) = N. The cluster-size distribution is defined 
as n(s)/ n(s). For J < J c , the cluster-size distribu- 
tion decays exponentially for large s. However, it decays 
in the power law form @ at the critical point J = J c , 
and the associated exponent r depends on the parame- 
ters 77 and A. We determine r using the generating func- 
tion method in the next subsection. For J > J c , a giant 
cluster forms and the distribution of finite-size clusters 
decays exponentially. The cluster-size distributions for 
various values of J are shown in Fig. ^ 



B. Generating function of the cluster-size 
distribution 



The probability that a vertex belongs to a cluster with 
size s is given by sn(s)/N, which is denoted as pis). In- 
voking the percolation theory, p(s) follows a power law 
with an exponential cutoff, 



p(s) 



-S I S c 



(38) 



where s c is the characteristic size, which depends on J 
and system size N. In the thermodynamic limit N — ► 00, 
s c diverges at J = J c . As in the percolation theory, 
the generating function V[z) = J2 s p( s ) zS 1S usem l f° r 
studying structural feature of the synchronized clusters, 
since its singular behavior is related to the critical be- 
havior of the synchronization transition, (i) The order 
parameter r ~ S/N can be obtained from the relation 
r ~ limjv^ooll - V(z* N )\, where P{z* N ) = E s <s-P( s ); 
i.e., the contribution by finite-size clusters. This can be 
achieved by choosing z* N ~ e _1 / Sm , where S m is a cluster 
size smaller than the largest cluster but larger the second 
largest cluster, (ii) From Eq. one can find that V(z) 
diverges for z > z c = lim s _»oo p(s)~ 1 ^ s , i.e., z c w e 1 / Sc . 
Thus, at J = J c , V(z) ~ (1 — z) T_1 as z — > z c = 1 in 
the thermodynamic limit. Thus, finding the singularity 
of P(z) enables one to obtain p(s). 

For the purpose, we introduce another generating func- 
tion Viz) as a partner of the local field f. From V(z), 
one can define a probability p(s) via the relation P{z) = 
^2 s p(s)z s , where p(s) is defined similarly to p(s) as the 
probability that a vertex belongs to a synchronized clus- 
ter of size s composed of the vertex and s — 1 neigh- 
boring vertices. For finite N, the generating function 
P(z) is analytic for \z\ < 1 and so is its inverse func- 
tion T~ l {z). To investigate the singularity of V(z) near 
z = 1, we consider the expansion of the inverse function 
z = V^iuj) = 1 - Ys n >\ M 1 - w) n around u = 1. The 
coefficient b n depends on J. Using Eqs. l(T3|) and 
and replacing f by 1 — u>, we can find that the generating 
function Viz) satisfies the self-consistent relations given 
below: 



z=viz)+Y^B n [j{i-nm 2n+1 

n=0 

+ C oc [J(l-P(z))]( A - 2 )/ f ' + .... 



(39) 



for J(l - P(z))k? n > 1, and 

OO 

z = Piz) + Y,A n [Jil-Piz))} 2n+ \ (40) 

n=0 

for J(l —Piz))^ <§C 1. Similarly, V(z) is determined as 

z-V{z) = £fl„[J(l-7>(z))] 2 " +1 

n=0 

+ C 00 [J(l-P(*))p- 1 >/' ? + ... 1 (41) 
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FIG. 1: (Color online) Cluster-size distributions n(s) / n(s) of synchronized clusters for the networks generated with A = 4.0, 
77 = 0.0, and N = 3000 at J = 0.7 J c (O), J = 0.8 J c (□)" J = 0.9 J c (A) for (a), J = J c (O) for (b), J = 1.1 J c (O), J = 1.2 J c 
(□), and J = 1.3J C (A) for (c), respectively. Solid lines drawn for reference have a slope of 7/3 for all. 



for J(l > 1, and 

OO 

z-7>(z) = 5>„[J(l-P(z))] 2 "+\ (42) 

ri=0 

for J(l--P(z))fc^ < 1. 

C. Behavior of p(s) at the critical point 

Here, we calculate the probability p(s) to find a vertex 
in the s-size cluster at the critical point J — J c explicitly 
in each case defined in SeclIIII 



In the case (I), since r\ < (A — 2)/3, we obtain that 

z = 1 + B\ (J c (l — 7^( z ))) 3 + ■•• to the leading orders by 
expanding z = V^iyj) around u> = 1 in cither Eqs. i|39|) 
or Thus, we obtain that 1 - V(z) - (1 - z) 1 / 3 , 

leading to 



p(s) 



-4/3 



(43) 



for large s. 

Using the obtained leading behaviors of V(z) around 
z = 1 in Eqs. (|41|) and Q42[l. one obtains the behavior of 
P(z) around z = 1 as 1 - "P(z) ~ (1 - z) 1 / 3 . Thus, 



p(s) 



-4/3 



(44) 



In the case (II), the singular term in Eq. I|39|) is 

relevant. In this case, the behavior of V(z) for z < z c 
differs from that for z > z c . z c is determined by the 
criterion J(l — V(z c ))k''^ n ~ 1. This case also happens for 
the cases (III) and (IV). 

In this case, the singular term C 00 [J c {l-V{z)))]^ 2 ^' 1 
is dominant in Eq. I|39|l : therefore, it follows that 1 — 
Viz) ~ J-\\ - zyil^-i) a t J = J c; which is valid for 
z ^> z c . From this result, p(s) is obtained as 



p(s) 



j-l s -(»7+A-2)/(A-2) 



(45) 



which is valid for s -C s c . 

On the other hand, when J(l — P(z))fc^ <§; 1 so that 
|A„ + i/A n |[J(l - P(z))] 2 < 1, one can obtain that 1 - 
P(z) ~ J- 1 |Ai|- 1 / 3 (l - z) 1 / 3 for z < z c ; therefore, 



(46) 



for large s » s c . s c is evaluated as follows: Substituting 
the result of 1— V(z) in the criterion J c (l — "P(z c ))fc^ ~ 1 
and using (1 — z c ) ~ s^ 1 , one can obtain system-size 
dependence of the characteristic size s c explicitly as 

Sc-jfc^-a^JV^-aW-D, (47) 

which diverges as iV — > oo. 

Together with Eqs. (|43|l and ijJSJl, we obtain that 



p(s) 



• s -(r )+ A-2)/(A-2) ( s « Sc ) j 
fc^- 2 )/ 3 -" S -V3 ( s » Sc ). 



(48) 



Next, using the result of 1 — V{z) ~ 1 — P(z) obtained 
from both Eqs. (|41|l and l|42|) . one can find that p(s) 
behaves similarly to p(s). That is, 

J s - (f7+ A-2)/(A-2) (s<<Sc)j 

p(s) ~ \^- a ^--v» ( s » Sc ). (49) 

In the case (III), the critical point J c is fi- 
nite in finite-size systems as being of order J c ~ 
k^ 2 -*! ~ iv(A-2-r,)/(A-i) pi U ggi n g the TV-dependence 
into Eq. I|4bj) and the expression, 1 — P(z) ~ 1 — z + 
C^lJcil - z)]( A " 2) /" for s < s c from Eq. one ob- 

tains p(s) as follows: 



p(s) 



• fc (A-2-,)(A-2)/, s _ (A _ 2+r))/l) (s <<; Sc); 



, -2(A-2)/3 4/3 
rum. o 



is > s c ). 



(50) 



Next, we derive p(s). We find that the leading singular 
term in Viz) for the case 1 — z s^ 1 ~ fc 2 n ~ A shows up 
in two ways. Substituting 1 — Viz) f» 1 — z + Coo[J c (l — 
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2)p-2)/i7jto Eq. gTJ), we obtain that 1-V(z) w 5 J C (1- 
^)+SoJ c C 00 [J c (l-z)]^- 2 )/"+C 00 [J c (l-z)]( A - 1 )/" + . • • . 
We compare the second with the third terms in order of 
magnitude. Using the fact that J c ~ k^ 2 ~ v , we find 
that there exist two subcases for s -C s c . The second 
term BoJ c Coc{Jc(l — z)]( x ~ 2 >'' n is more dominant than 
the third term C^J^l - z)]^" 1 )/'' when 1 - z < s" 1 
and vice versa. Here, it is found that a new crossover size 
s* scales as 

s.~fc<r i ><"- A+a >. (si) 

From the behaviors of T'(.z) in the three different sub- 
cases, we obtain the probability p(s) as 

p(«) ~ < fc(^ 2 -"» A - 2+r ')/" s -(^- 2 +")/" ( s , « a « Sc ), 
lfc^ 2)/3 -" S - 4 / 3 (s»s c ). 

(52) 

One can notice that the subcase s -C s* diminishes when 
7? < 1, but it is extended as the parameter ry increases. 

In the case (IV), the third term C 00 [J c (l-z)]( A - 1 )/ ? ? 
in 1 — V{z) in the case (III) is always dominant when 
l-z> s^ 1 . Moreover, A in Eq. ij^l diverges as A ~ 
^wT j which has to be considered in the relation, 1 — 
P(z) ~ AjJ c (l - for 1 - 2 < s™ 1 . Consequently, 

p(s) behaves as 




-1 1/3 1 2 



il 

FIG. 2: (Color online) Diagram in the space of (rj, A) of eight 
different domains each corresponding to a distinct synchro- 
nization transition. The transition nature of each domain is 
listed in Table Q 

V. LARGEST CLUSTER SIZE AND 
FINITE-SIZE SCALING 

In this section, we first investigate the iV-depcndence 
of the largest cluster size at J c . Next, based on this 
result, we derive a finite-size scaling form for the order 
parameter near J c . 



p{s) 



• fc (A-l)(A-2-„)/„ s _ (A _ 1+r))/r) (s <<; Sc)) 



Jl-3A)/3 s _4 /3 



(s > S c ). 



(53) 



To substantiate the predictions of this section, we in- 
vestigate the asymptotic behavior of p(s) in a numerical 
manner. The static model introduced in |llj| is used for 
underlying network in our simulations. The network has 
TV = 3000 oscillators and its mean degree (k) is 4.0. The 
values of A and i] are chosen as 4.0 and 0.0, respectively. 
This pair belongs to the case (I). First, we simulate the 
system at J = J c to determine Cth defined in Sec. II V A1 
During the simulation, we assume a large value of Cth 
and then collect the pairs of vertices that the Cy of each 
pair is larger than the assumed Cth- After that, we de- 
termine clusters and obtain the cluster-size distribution. 
We then adjust Cth by somewhat decreasing or increas- 
ing it, and repeat these procedures until the power-law 
distribution appears in the cluster-size distribution. If 
the cluster-size distribution follows the power-law form 
of p(s) ~ s _r , the corresponding value of Cth is consid- 
ered as the threshold value Cth- It is found numerically 
that Cth ~ 0.7, independent of the system size N. In 
our simulations, we obtain r + 1 ks 7/3, which is close to 
the theoretical value in Eq. (|44|l . as shown in Fig. ^b). 
We also performed simulations for various J < J c and 
> J c as shown in Figs. Ufa) and (c), respectively. The 
power-law behaviors in the cluster-size distribution do 
not appear in these cases. 



A. The largest cluster size 

The largest cluster size S can be obtain from the rela- 
tion [m, 



s>S 



S_ 
N' 



(54) 



In the case (I), we use the result of Eq. JUJl, and 

obtain simply that 



S - N 3/i . 



(55) 



In the case (II), p(s) displays a crossover at s c , and 
thus the obtained value of the largest cluster size must 
satisfy the self-consistency conditions. For instance, the 
largest cluster size obtained by Eq. (|54|) for s -C s c in 
Eq. H49|) has to be smaller than s c . As a result, the largest 
cluster size behaves differently in the two subcases r) < 1 
and r\ > 1, which we denote (Ua) and (lib), respectively 
In each subcase, we obtain that 



S 



AT(4A-5-3„)/[4(A-l)] in (Ha)) 
N (X-2)/(X-2 +n ) in (nb ). 



(56) 



The largest cluster size S in (Ha) was determined from 
p(s) for s 3> s c , and is indeed much larger than s c , 
whereas it in (lib) was done from p(s) for s <C s c . 
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TABLE I: The probability to find a vertex in s-size cluster p(s), the largest cluster size 5 at the critical point, and the critical 
exponents (3 and fx for the eight cases shown in Fig. 



domain 


p(s) 


s 


/3 


M 


(I) 


s -4/3 


AT3/4 


1 

2 


2 


(Ha) 


/ s-<*+»-W-*) ( a « Sc ), 


jy(4A-5-37,)/[4(A-l)] 


?7 


V 4(A-1) 
1 + 3?7 A - 2 - r/ 


(lib) 




Ar (A-2)/(A-2+^) 


X 9 ti 


A - 2 + J? 
A - 2 - 77 


(Illal) 




k (X-2- v)i X-2 +v) / Vs - (X -2 +v)/v ^ <<; s <<; 
fc (A- 2 )/3^ s -4/3 (g ^ 


A f(4A-5-3^)/[4(A-l)] 




7? 4(A-1) 
1 + 377 77 - A + 2 


(IIIa2) 


1 


V 




(nib) 




jy J 7/(A-2+7 7 )+(A-2-»7)/(A-l) 


V - A + 2 


277 

77- A + 2 


(IIIc) 




Ar (A-2)/(A-l+^) 






(IV) 




A T(A-2)/(A-l+^) 


A- 1 


2(A-1) 


\k^ 3 s-^ (s» So ). 


r?- A + 2 


77- A + 2 



In the case (III), p(s) exhibits three distinct power- 
law behaviors. Thus, this case is divided into three sub- 
cases. They are as follows: 77 < 1 (Ilia), 1 < 77 < 
\/A 2 - 3A + 3 (Hlb), and rj > VA 2 - 3A + 3 (IIIc). The 
largest cluster size in each subcase is given as 



' iV (4A-5-3r ; )/[4(A-l)] 
{ N r,/(^~2+ n ) + (X~2~ V )/(X-l) 

N (x-2)/(\-i+ n ) 



in (Ilia), 

in (Illb), (57) 

in (IIIc). 



In the case (IV), the largest cluster size is deter- 
mined simply by p(s) for s <C s c since the resulting 
largest cluster size fulfils the criterion S < s c for r\ > 0. 
Thus, 



5- ^ jv( A -2)/(A-l+r,)^ 



(58) 



for the cases (I) and (II), and 

r = N-P^iPiJN 1 ^) 
for the cases (III) and (IV), where 



const for iCl, 



for x ^> 1. 



,/3 



(60) 



(61) 



The critical exponent \x is determined by the relation 
r c - 7V-/Vp. The value of [i varies depending on the 
cases determined by the magnitude of 77 and A. 

We present the diagram in Fig. [21 comprising eight dis- 
tinct cases in the (77, A) plane. Each case in the dia- 
gram corresponds to a distinct behavior of the critical 
exponents /3 and /j, the cluster-size distribution, and the 
largest cluster size. We summarize those features in Ta- 
ble □ 



B. Finite-size scaling 

Here, we evaluate the magnitude of the order parame- 
ter r c at J c and establish the finite-size scaling function. 
To proceed, we compare the magnitude of cooperative 
synchrony S/N with that of the background synchrony 
~ N^ 1 / 2 . The order parameter r c is defined as ~ S/N if 
S/N > iV~ 1/2 , and ~ N~ 1/2 otherwise. Under this crite- 
rion, we obtain r c as ~ S/N in the cases (I) and (II), and 
~ 7V-V2 in tnc cases (Illb), (IIIc), and (IV). The case 
(Ilia) is divided into two subcases, 2A — 3i] — 3 > (<)0. 
They are denoted as (Illal) and (IIIa2), respectively. 
The order parameter r c behaves as ~ S/N and ~ N^ 1 / 2 
in (Illal) and (IIIa2), respectively. 

By using that r ~ A^ 3 and A^-dependent behavior of r c 
at J c , we can construct a finite-size scaling form as 

r = N'^^^AN 1 ^) (59) 



VI. NUMERICAL SIMULATIONS 

We perform direct numerical integration of Eq. J2J to 
confirm the analytic results. In particular, the finite-size 
scaling behaviors in Eqs. (|59|l and l|6U|) are compared. For 
the purpose, we generate random SF networks using the 
static model [ill witn system sizes of A r = 400, 800, 1600, 
and 3200, mean degree of (k) = 4.0 and values of 77 and 
A chosen from each domain in Fig. [5J Numerical values 
of (A, 77) we used are listed in Table ITT1 For the numerical 
integration, we apply the Heun's method fl2T ]. Time is 
discretized by a unit step St = 0.005 and runs up to 
t = 1.2 X 10 4 . Ensemble average is taken over O(10 2 ) ~ 
O(10 3 ) different configurations of natural frequencies and 
network realizations, respectively. 

Numerical results are presented in Fig. |3| For each 
Fig. [3f a)— (c) , the critical point J c is finite. We find J c 
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TABLE II: Numerical values of the parameters (rj, A) we used for Figs. [3] /3t and \it are theoretical values for a given set of (77, A) 
in the third column. f3 n and are numerical values to draw Figs for each case. For (a)-(c), the theoretical and numerical 
values are the same each other for both /3 and fj,. However, they can be different for (d)-(h). rt and r n are the order parameters 
in scaling form formulated with the theoretical values of (3t and [i t , and the numerical values (3 n and fi n , respectively. 



Fig.! 


domain 


(V, A) 


A 


0n 




fJ-n 


r t 


r„ 


(a) 


(I) 


(1/3,4) 


1/2 


1/2 


2 


2 


N~ 


-V^CAiV 1 / 2 ) 


N- 


"^(AA^ 2 ) 


(b) 


(Ila) 


(5/6,4) 


5/7 


5/7 


49/120 


49/120 


AT 


" 7 / 24 V(AiV 49 / 120 ) 


N~ 


- 7 / 2 V(AA^ 49/120 ) 


(c) 


(lib) 


(4/3,4) 


2 


2 


5 


5 


N~ 


' 2 / B ip(AN 1/& ) 


N- 


- 2 /V(AAr 1/5 ) 


(d) 


(Illal) 


(1/3, 13/6) 


2 


20/7 


14/3 


20/3 


N~ 




N- 


- 3 /V(JiV 3 / 20 ) 


(e) 


(IIIa2) 


(2/3, 13/6) 


4/3 


50/23 


8/3 


100/23 


N~ 




N- 


"^^(JiV 23 / 100 ) 


(f) 


(nib) 


(5/2,4) 


5 


5 


10 


10 


N~ 


- 1 / 2 V(Jiv 1/10 ) 


N~ 


"^^(JiV 1 / 10 ) 


(g) 


(IIIc) 


(11/4,4) 


11/3 


11/3 


22/3 


22/3 


N~ 


- 1 ' 2 4,(JN 3 ' 22 ) 


N- 


-^XiA^ 22 ) 


GO 


(IV) 


(4,4) 


3/2 


5 


3 


10 


N~ 


- 1/2 i>(JN 1/3 ) 


N' 


-^(JA^ 10 ) 




(b) 


1 1 GJK 

i 




(Ila) 



-5 5 

. ,,49/120 



. (f) 


1 

□ 
□ 

6 • 

□_o 

[% 
CO 




(1Mb) " 






■ 1 r 

□ L 


" (g) 


O 

□ o - 








O A 




u* 








□ 




.00 




(NIC) - 







6 
5 
4 

\ 3 
2 

1 





1 1 

- (d) 


1 1 1 a t 

* 

"o 
ao 
□0 

* 






1 1 


<r 
<P 

§ (Illal) - 



0.5 1 1.5 



<3 




0.02 0.04 



0.02 0.04 0.06 

jn" 22 



0.001 0.002 0.003 0.004 0.005 



FIG. 3: (Color online) Finite-size scaling behaviors of the order parameter r. Data are collected from the static network model 
with mean degree (k) = 4 and system sizes N =400(O), 800(A) 1600(O), an d 3200(D). Numerical values of the tunable 
parameters (77, A) and the critical exponents are given in 

Table. E] for each case. For the critical point J c , theoretical ji*-* and 
numerical J< n) values used in (a)-(c) are different as (J c (t) , ji n) ) = (0.92, 1.32) (a), (0.37,0.50) (b), and (0.13,0.18) (c). 



numerically to make the obtained numerical data col- 
lapsed for different system size N in the scaling plot with 
theoretical values of (5 and fi. Theoretical and numeri- 
cally found values of J c , denoted as and ji n \ respec- 
tively, are compared as (Jc , Jc ) = (0.92, 1.32) for (a), 
(0.37,0.50) for (b), and (0.13,0.18) for (c). They belong 
to the cases (I), (Ila), and (lib), respectively. 

For Figs. 3 (d)-(h), the critical point J c is zero. For 
the cases of (d),(e), and (h), we find that numerical data 



do not collapse well in the scaling plot of rN^^ versus 
JN 1 ^ with theoretical values of fit and /it tabulated in 
Table I. Instead, we adjust numerical values of /?„ and fi n 
values empirically to make the obtained numerical data 
collapsed. Those empirical values of (3 n and fx n are com- 
pared with the theoretical values as listed in Table II. 
The discrepancy in (d) and (e) originates from the pres- 
ence of intrinsic degree-degree correlation in the static 
model when the degree exponent 2 < A < 3, while the 
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theoretical values were obtained under the assumption 
that the degree-degree correlation is absent. 



VII. CONCLUSIONS AND DISCUSSION 

In this paper, we have investigated the nature of the 
synchronization transition generated by -/V limit-cycle 
oscillators located at random SF networks with degree 
exponent A. The dynamics is given by a modified 
Kuramoto equation with the asymmetric and degree- 
dependent weighted coupling strength in the form of 
Jk\ *', where fcj is the degree of vertex i. Depending 
on the sign and magnitude of 77, the influence of the hub 
vertices on the dynamics can be moderated or amplified, 
and determines the nature of the synchronization transi- 
tion. Applying the mean-field approach to the modified 
Kuramoto equation, we derived the critical point, the 
size distribution of synchronized clusters, and the largest 
cluster size at the critical point. The critical exponents 
associated with the order parameter and the finite-size 
scaling are determined in terms of the two tunable pa- 
rameters, (77, A). All results are summarized in Table. [I] 
The parameter space of (77, A) is divided into eight dif- 
ferent domains, in each of which the transition nature is 
distinct. 

It would be interesting to notice that the critical ex- 
ponents (3 and (i associated with the order parameter 
and the finite-size scaling of the synchronization tran- 
sition depend on the parameter 77 defined in the cou- 
pling strength. The result is unusual from the perspec- 
tive of the universality in the critical phenomena in reg- 
ular lattice where the details of the cou plin gs are mostly 
irrelevant unless they are long rang ed hj. This result 
implies that structural features of SF networks such as 
the degree distribution are not sufficient to understand 
the dynamic process on such networks. Asymmetric cou- 
pling in dynamics is a relevant perturbation in such net- 
works, because of the heterogeneity of the degree distri- 
bution. Such a behavior was also observed in the sandpile 
model Ha. 
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APPENDIX A: DERIVATION OF EQS. (1l9l2l)l> 

In this Appendix, we present the derivation of 
Eqs. (|19H20|) from Eqs. O and US with 77 > 0. Since 
the right-hand-side of Eq. I|ltjfl becomes the same as that 
of Eq. (|15f) when A is shifted by 1, we here analyze in 



detail the nature of Eq. I|15fl only, which applies also to 
Eq. l(Tfi|) with A replaced by A + 1 . 

The coefficients in Eq. I|15fl diverge with increasing N 
for n > n c = |~(A — 77 — 2)/ (277)] , where \x] is the smallest 
integer not smaller than x, due to the divergence of the 
generalized harmonic number for < q < 1. While 
H^ n ~ £(g) for q > 1, H^ n diverges in the limit m — > 00 for 
< q < 1 and its asymptotic expansion can be obtained 
by using the relation = ((q) — ((q,m + 1) and the 
asymptotic expansion of the Hurwitz zeta function [T^j 



C(<?,rn) 



1 



,1-9 



dx 



sin((7 tan 1 x) 



(l + x 2 )i/ 2 (e 27Tmx - 1)' 



(Al) 



One can see that the integral is of order m 1 in the limit 
m — > 00 since 

sin(<7 tan -1 x) 



dx - 

o (l + x 2 )i/ 2 (e 2 ™* -1) 

l/m 

< / dx— - — h / dx 



1 



a27tmx 1 

e 1 J l/m 

0(m _2 )+0(ro _x ) J 
and thus for < q < 1, behaves as 
(m + 1) 1 ~« 



(A2) 



k=l 



k'i 



1-9 



+ C(«z)+0(m-«). (A3) 



When 77 > 0, the terms with such diverging coefficients 
exist and thus we can rearrange the expansion as follows: 

~ ( n - l/2)l(-l) n (Jf) 2n+1 H*- v - 2nv - 1 



n=0 



2»+8/a n !(n + l)!C(A-l) 



(A4) 



= f> (rc-l/2)!(-l) n C(A- V -27777.-1) 2n+1 
2^ 2«+ 3 / 2 77!(n + l)!C(A-l) 1 ' 

n=0 • ' x ' 

« ( n _ 1 /2)!(-l)"fc 2 +"+ 2 "'7- A (Jf) 2 "+ 1 
+ ^ 2™+ 3 / 2 n!(n + 1)!(2 + 7; + 27777 - A)C(A - 1) 

OO 

= Y, B n (Jf) 2n+1 + (Jf)^^C(Jfkl), (A5) 

n=0 

where we approximated H^ _1 by £(A — 1) since A — 1 > 
1. The coefficients B n are defined in Eq. I|17|) and the 
function C(x) is defined by 

C(x) = c n x {2+ ' n+2mi - x)/ \ with 

n—n c 

(n- l/2)!(-l)« 

Cn ~ 2 n +*/ 2 n\{n + 1)!(2 + 77 + 2n?7 - A)C(A - 1) ' 

(A6) 

While C(x) behaves as x ( - 2+r i- x + 2 ^)/n f or x < 1, it 
converges to a constant Coo for x — > 00 yielding a non- 
analytic term C 00 ( Jf)( A_2 )/ r ' as can be seen in Eq. Q22p. 
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Therefore, the magnitude of Jfk^ is essential for the 
determination of the leading behaviors of the right-hand- 
sides of Eqs. (fTo)> and l(TH|l for_f <C 1. If Jfk 1 ^ > 1, one 
can approximate the function C( Jffc^) by a constant Coo 
that is evaluated as 



lim 



oo 



(n ~ l/2)\(-l) n 2-<- x+2r i-V/( 2 n) 
x^oo^j n \(n + 1)!(2 + 77 + 2nr? - A)C(A - 1) 

X (2- 1 /2 a ,)(2+>)+2n»,-A)/7 ) 

1 ^ (-l) n (n- 1/2)! 

" 2( A + 2 ')-2)/(2>))^(^ _ !) 

^ ' n—n c 



1+77+2™^- A 



n!(n + 1)! 



< j dyy 1 

[(A- 2^- 2)7277]! [(2- A-t?)^]! 
7?2 (A+4^-2)/2^( A + ^ _ 2 ) /27y]!C(A - 1) ' 



(A7) 



The function C(x) defined in the text can be approxi- 
mated in the same way by a constant that is identical 
to Coo with A + 1 in place of A. Therefore, one should 
refer to Eqs. 1-21) and (|23[l for the correct expansions of 
r and f around f = in the case of Jfk 7 ^ 3> 1 while 
Eqs. (|15(l and Q16|l can be used in the case of Jrk^ n <C 1. 
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